Week 11: Linear Modeling
When analyzing a single numerical variable, we can use numerical statistics to gain insights into key attributes such as
Alternatively, we can use graphical tools like box plots and density plots to represent these characteristics.
maxt <- read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%
pull(maxt)
p1 <- ggplot() +
geom_boxplot(aes(maxt)) +
xlab(expression(degree*C)) +
theme_light()
p2 <- ggplot() +
geom_density(aes(maxt)) +
xlab(expression(degree*C)) +
theme_light()
patchwork::wrap_plots(p1, p2) +
patchwork::plot_annotation(title = "Maximum daily temperature for bushfires in Victoria between 2015 and 2018")When working with two numerical variables, a scatterplots is a natural choice to explore association between them.
It helps identify whether the variables are:
read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%
ggplot() +
geom_point(aes(mint, maxt), alpha = 0.6) +
xlab(expression(Minimum~temperature~(degree*C))) +
ylab(expression(Maximum~temperature~(degree*C))) +
theme_light() +
ggtitle("Daily maximum and minimum temperature for bushfires \nin Victoria between 2015 and 2018")Correlation is the linear association between two variables, ranging from \(-1\) to \(1\).
For two variables, \(X\) and \(Y\), the correlation coefficient \(r\) is calculated as follows:
\[r_{xy} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n}(y_i - \bar{y})^2}} = \frac{S_{xy}}{S_x S_y},\]
where \(S_{xy}\) represents the biased sample covariance between \(X\) and \(Y\), and \(S_x\) and \(S_y\) denote the biased standard deviations of \(X\) and \(Y\), respectively.
Tip
Interpretation:
Causation means that one variable directly influences another.
Here are some examples:
Remember: Correlation does not imply causation!
Modeling is a method for understanding the relationship between different variables.
We typically represent the relationship of one response variable, \(Y\), in relation to other explanatory variables, \(X_1, X_2, \ldots, X_p\), using a mathematical function:
\[Y = f(X_1, X_2, ..., X_p).\]
Why we need modeling?
While we can observe the data, we may not have insights into the underlying processes, i.e. the joint distributions of the variables.
Therefore, we rely on data to make inferences about these processes and build models to explain them.
There are lots of different models.
For now we focus on linear models.
read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%
ggplot() +
geom_point(aes(mint, maxt)) +
geom_smooth(aes(mint, maxt), method = "lm", se = FALSE) +
xlab(expression(Minimum~temperature~(degree*C))) +
ylab(expression(Maximum~temperature~(degree*C))) +
theme_light() +
ggtitle("Daily maximum and minimum temperature for bushfires \nin Victoria between 2015 and 2018")A simple linear regression (SLR) is a regression model that involves a single explanatory variable. The formula is expressed as follows:
\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad i = 1, \ldots, n.\]
ggplotWhat are the intercept and slope?
ggplotWhat are the intercept and slope?
lm\[\text{maxt}_i = \beta_0 + \beta_1\text{mint}_i + \varepsilon_i, \quad i = 1, \ldots, n.\]
lm\[\text{speed}_i = \beta_0 + \beta_1\text{dist}_i + \varepsilon_i, \quad i = 1, \ldots, n.\]
But how are these coefficients actually estimated?
First, we need to introduce the concepts of fitted values and residuals.
\[\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i.\]
The \(i\)-th residual \(e_i\) is the difference between the fitted value and the actual value of the response variable:
\[ e_i = y_i - \hat{y}_i = y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i). \]
We aim to minimize the residuals, meaning we want the regression line to be as close to the data points as possible.
But how do we determine which solution is better? How do we evaluate it?
We can minimize the mean squared error (MSE) by solving for the coefficients:
\[ \underset{\beta_0, \beta_1}{\text{argmin}} \, \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \]
This ensures that the regression line would be the line closest to the cloud of points!
The minimization process described before leads to the following estimates:
\[\hat{\beta}_1 = \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n(x_i -\bar{x})^2} \Rightarrow \text{slope}, \text{ and}\]
\[\hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x} \Rightarrow \text{intercept},\] where \(\bar{y} = \frac{1}{n}\sum_{i=1}^ny_i\) and \(\bar{x} = \frac{1}{n}\sum_{i=1}^nx_i\) are the sample mean.
Note
Recall that \[r_{xy} = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n}(y_i - \bar{y})^2}}.\]
So, for simple linear regression, \(\hat{\beta}_1 = r_{xy}\frac{S_y}{S_x}\).
Once the model is fitted and all parameters are estimated, we can make predictions by plugging in \(x\) to obtain the predicted value \(\hat{y}\).
This predicted value reflects the average \(y\) for each corresponding value of \(x\), i.e. \(\hat{E}[y|x]\).
\[\widehat{\text{maxt}}_i = 12.02 + 1.04 \times \text{mint}_i.\]
Examples:
Warning: While we “expect” this to happen, there will be some variability, which is captured by the error term.
Note that, when \(\text{mint} = 0\), the predicted \(\widehat{\text{maxt}} = 0 \times 1.04 + 12.02 = 12.02\).
Similarly, the slope \(\hat{\beta}_1 = 1.04\) can be interpreted as follows:
We have discussed the most basic regression model in statistics. Now, let’s consider its extension to include multiple explanatory variables:
\[y_i = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... +\beta_px_p + \varepsilon_i, \quad i=1,...n,\]
where \(\varepsilon_i \sim N(0, \sigma^2)\) for \(i = 1,...n\).
Model assumptions:
In a CNLRM, the interpretation of the intercept term \(\hat{\beta}_0\) remains the same.
For the other parameters \(\hat{\beta}_j\) where \(j = 1, \ldots, p\), the interpretation is as follows:
All models are wrong; some are useful - George Box.
After fitting a linear model, it is crucial to assess whether the model assumptions are satisfied.
There is a real risk, that a model is imposing structure that is not really there.
It can result in inaccurate estimates of the coefficients and misleading interpretations of the relationship between \(x\) and \(y\).
The first thing to check is the residuals versus fitted values plot.
What constitutes a good residual plot?
The second thing to check is the residuals versus each individual explanatory variable plot.
The principles for a good residual plot remain the same.
This type of plot offers clearer insights into the relationship between the residuals and each individual explanatory variable, which can be invaluable for refining and improving the model.
The third thing to check is the Q-Q plot.
To check a Q-Q plot, compare the distribution of your residuals to a theoretical normal distribution.
In a good Q-Q plot, the points should closely follow the 45-degree reference line.
Significant deviations from this line, such as curves or outliers, indicate that the residuals may not be normally distributed, which could suggest violations of model assumptions.
We can display the actual residual plot alongside a matrix of residual plots simulated under the assumption that the model is correct.
If the true residual plot is distinct and identifiable within the matrix, it provides evidence suggesting a violation of the model assumptions.
This method is known as the lineup protocol.
Beyond checking model assumptions, we often need to select between different models. For instance, choosing between maxt ~ rf + mint or maxt ~ mint.
There are several statistics to assess model fit:
\[R^2 = 1-\frac{\sum_{i=1}^{n}(y_i-\hat{y})^2}{\sum_{i=1}^{n}(y_i-\bar{y})^2} = \frac{\text{Residual sum of square}}{\text{Total sum of square}} \in [0,1].\]
Which model is better?
AIC, BIC, and Deviance are goodness of fit measure to compare models.
AIC = Akaike Information Criterion (can be used to compare models. The smaller the value the better the model.)
Similarly BIC = Bayes Information Criterion indicates how well the model fits, best used to compare two models. Lower is better.
Deviance is the residual variation, how much variation in response that IS NOT explained by the model. The close to 0 the better, but it is not on a standard scale. In comparing two models if one has substantially lower deviance, then it is a better model.
Hans Rosling was a Swedish doctor, academic and statistician, Professor of International Health at Karolinska Institute. Sadly he passed away in 2017.
He developed a keen interest in health and wealth across the globe, and the relationship with other factors like agriculture, education, energy.
You can play with the gapminder data using animations at https://www.gapminder.org/tools/.
gapminderContains subset of the data on five year intervals from 1952 to 2007.
nestETC1010/ETC5510 Lecture 11 | Melbourne time